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We find numerical solutions of the coupled system of Einstein-Maxwell's equations with a linear 
approach, in which the magnetic field acts as a perturbation of a spherical neutron star. In our 
study, magnetic fields having both poloidal and toroidal components are considered, and higher order 
multipoles are also included. We evaluate the deformations induced by different field configurations, 
paying special attention to those for which the star has a prolate shape. We also explore the 
dependence of the stellar deformation on the particular choice of the equation of state and on the 
mass of the star. Our results show that, for neutron stars with mass M = 1.4 Mq and surface 
magnetic fields of the order of 10 15 G, a quadrupole ellipticity of the order of 10~ 6 — 10~ J should 
be expected. Low mass neutron stars are in principle subject to larger deformations (quadrupole 
ellipticities up to 10~ 3 in the most extreme case). The effect of quadrupolar magnetic fields is 
comparable to that of dipolar components. A magnetic field permeating the whole star is normally 
needed to obtain negative quadrupole ellipticities, while fields confined to the crust typically produce 
positive quadrupole ellipticities. 

I. INTRODUCTION 

The measured periods and spin down rates of soft-gamma repeaters (SGR) and of anomalous X-ray pulsars (AXP), 
and the observed X-rayluminosities of AXP, indicate that these neutron stars have extremely high magnetic fields, as 
large as 10 14 — 10 15 G [J 0] . Furthermore, if these sources are the central engine of gamma-ray bursts, as suggested in 
[H, their magnetic field might even be larger. Up to now, about ten highly magnetized neutron stars, the "magnetars" , 
have been identified in our Galaxy, but their actual number may be larger, and it has been suggested that a fraction 
of pulsars ( > 10% would possibly become magnetars at some stage of evolution. The discovery of magnetars has 
triggered a growing interest in the study of the structure, dynamics and evolution of neutron stars with large magnetic 
fields, and has raised a number of interesting issues. For example, quasi-periodic oscillations have been detected in 
the aftermath of the giant flares of SGR 1806-20 and SGR 1900+14, and it is not clear whether they are associated 
to crustal modes, or to modes of the magnetic field (or both); if the spacing between the observed frequencies would 
be explained, one may gain information on the internal structure of the star 

In addition, magnetars may be interesting sources of gravitational waves, especially if they possess a toroidal 
magnetic field; indeed, as suggested by Jones and Cutler Jy, 0], a large toroidal component tends to distort the star 
into a prolate shape, leading to a secularly unstable object: the wobble angle between the angular momentum and 
the star's magnetic axis would grow on a dissipation timescale, until they become orthogonal. This may produce a 
copious flux of gravitational waves, potentially detectable by the advanced version of gravitational wave detectors 
LIGO and VIRGO @]. 

In order to understand magnetars' structure and dynamics, it is necessary to model their equilibrium configuration 
in the framework of general relativity, including both poloidal and toroidal magnetic field components. Toroidal fields 
should form during the first seconds after core collapse, when the star is likely to be rapidly and differentially rotating: 
the fluid motion would drag the poloidal field lines creating large toroidal fields Q; in addition, convective motions 
prevailing in the early life of a neutron star could also create toroidal fields by dynamo processes 0, H, Q. These 
toroidal components are expected to survive when the proto-neutron star cools down and the crust forms. We also 
stress that large toroidal components contribute to explain the giant flares in current models of SGR's [l[. In [l(J it 
has been shown that a purely poloidal magnetic field is unstable, and decays on a timescale much shorter than the 
star's life (see also [ll[ and references therein, and [13); however, as discussed in [TTj] , a magnetic field configuration 
with prevailing toroidal component is also expected to be unstable on a short timescale. Thus, both toroidal and 
poloidal magnetic fields have to be included to construct accurate, and stable, models of magnetars. 

In recent literature, magnetars equilibrium configurations have been studied by solving Einstein-Maxwell's equa- 
tions, coupled with the Hydrodynamics equations, in full general relativity [l3l [Til [l5| . However, the numerical 
schemes used in most cases require circularity of the space-time, i.e. the existence of two hypersurface-orthogonal 
Killing vectors, and this assumption automatically excludes toroidal magnetic fields, since they break circularity. 
Therefore, in [ID, [lj|, [l~5| only poloidal magnetic fields have been considered. 

A different approach has been used in [ilil [TtI ITsj . where equilibrium configurations have been studied using 
a perturbative techniques, i.e. solving Einstein-Maxwell-Hydrodynamics equations, linearized about a spherically 
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symmetric background, and expanding the perturbed equations in tensor harmonics. Toroidal fields have been included 
in the analysis only in [16[, but this work is based on very restrictive assumptions: the magnetic field is assumed to 
vanish outside the star. Poloidal and toroidal fields have also been considered in the framework of Newtonian gravity 
in a recent work 

In this paper we construct equilibrium configurations of neutron stars with strong magnetic fields, in general 
relativity. Since magnetars rotate very slowly, we restrict to non-rotating stars. However, rotation can play an 
important role in the early phases of the stellar evolution, therefore it will be included in future developments of this 
work. We follow a perturbative approach, generalizing the work of [ItJ to include toroidal magnetic fields, with a 
magnitude comparable with that of the poloidal fields. We start solving the relativistic Grad-Shafranov equation, to 
which Maxwell's equations can be reduced, in the background of a non rotating star; we impose a set of boundary 
conditions which correspond to different magnetic field configurations, and construct the corresponding stress-energy 
tensor. The magnetic field perturbs the star, which is consequently deformed; to compute the stellar structure and its 
deformation, we then solve the Einstein-Maxwell- Hydrodynamics equations linearized about the spherically symmetric 
background of the non rotating star, having the electromagnetic and the fluid stress-energy tensors as a source. We 
compare the deformation induced by a magnetic field with that which would be produced by rotation, and find that 
effect of magnetic fields is dominant for magnetars as SGR and AXP. We discuss how the magnetic field profile and 
the corresponding stellar deformation depend on the stellar mass and on the equation of state of the fluid composing 
the star, comparing different stellar models. In current literature, only the I = 1 multipole of the electromagnetic 
potential is usually considered. In this paper we also solve the relevant equations for the I — 2 multipole. 

The main features of the perturbative approach are described in Section[H] The results of the numerical integrations 
of the relativistic Grad-Shafranov equation, and of the equations of stellar perturbations, are reported and commented 
in Sect ion HTT1 for different field configurations and different stellar models. In Section Hill we also discuss the effects of 
the / = 2 multipole. Conclusions are drawn in Section HVl 



II. STRUCTURE OF A STATIONARY, AXISYMMETRIC NEUTRON STAR WITH POLOIDAL AND 

TOROIDAL MAGNETIC FIELDS 

In what follows we shall assume that the magnetized fluid composing the non rotating neutron star can be described 
within the framework of ideal magneto- hydrodynamics (MHD), i.e. that there is no separation of charge currents 
flowing through the star. It should be mentioned that, although this assumption is appropriate inside the fluid core, 
it may not apply to the stellar solid crust. The magnetic field, and the deformation it induces on the star, are treated 
as stationary and axisymmetric perturbations of a spherically symmetric background. We consider perturbation s up 
to order 0(B 2 ). We shall follow the notation and the formalism introduced by Konno, Obata and Kojima in pTj . 
generalized to include toroidal magnetic fields. 

Before proceeding with the perturbative approach, in the next subsection we shall summarize some general properties 
of stationary, axisymmetric magnetized stars, which will be useful in subsequent sections. These properties and their 
proofs can be found in the literature, but are scattered in different papers (20| . [2l| : here we report them in a unified 
and consistent way. 



A. Some properties of stationary, axially symmetric magnetized stars 

We consider a stationary, axisymmetric space-time describing a magnetized star, with coordinates 

= (t,x a ,c/)) (a =1,2) (1) 

where r\ — d/dt and £ = d/d(j> are Killing vectors. The coordinates x a can be, for instance, spherical coordinates 
(r, 6), or cylindrical coordinates (r, z). Any stationary, axisymmetric quantity, such as the vector potential or the fluid 
4- velocity, are independent of t and <f>, i.e.: = A^(x a ), — u^{x a ). 
The electric and magnetic field are defined as 

E M = F^u" , = ~e a ^ B u fi F^ s , (2) 

where F M „ = d^A v — d^A^ It may be noted that since d t A^ — d^A^ — 0, F t< p — 0. We define the local angular velocity 
of the fluid as 
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The components u a of the fluid velocity are called meridional currents. Furthermore, we define the quantities 

<S>(x a ) = rfA^ = A t , *(x a ) = £M„ = . (4) 
The ideal MHD hypothesis implies that the electric field, measured by a comoving observer, vanishes: 

E» = F^u v = . (5) 
An axially symmetric magnetic field is poloidal if its non vanishing space-components arc 

(B a ,0); (6) 

it is toroidal if 

(0,0,B*). (7) 

1. Vanishing meridional currents 
If meridional currents vanish u a = 0, and eq. gives 

E a = FatU 1 + FafU* = {F at + SlFarfu* =0 =► 
As F Qt = <9 Q <i> and F a = 9 \I>, eq. (JSJ) becomes 



Fit _ Fitj, 
F2t F24 



-CI. 



Furthermore, from (fST) we have 



d a <s> = -nd a * . 



di^* - d 2 $di* = , 



(8) 

(9) 
(10) 



which implies (assuming the domain where $ and VP are defined is simply connected) that $ = $(^). From eq. © 
then it follows 



d<I> 



-O. 



(11) 



As shown by Carter (see theorem 7 of [20} , and its corollary; see also [22|), if the space-time is stationary and 
axisymmetric, and if meridional currents are zero, then A a = 0. Therefore, the vector potential is 



^ = ($,0,0,*), 



the electromagnetic tensor becomes 



F, 


















o 


-*,a 


o 



(12) 



(13) 



and the magnetic field is 

B a = e a ^ u u F uu = (O,B a ,O) , (14) 
since g at = g a ^ = when u a = [20]. Thus, if meridional currents vanish the magnetic field is poloidal. 



2. N on-vanishing meridional currents 
Let us now consider the general case u a ^ 0. eq. © gives 



E t = F ta u a = -u a 9 a $ - -u^S = 



d<P 

dt> 



0, 



(15) 
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i.e., $, 'J are constant along the fluid flow. Then 

u 1 ^.! + u 2 $ 2 = 

U 1 *! +M 2 *, 2 = 0, (16) 

which implies $ = <&(^E'). We introduce the quantity 



so that 



^, (17) 



<9 a $ = -Q9 a * . (18) 

Notice that in general f2 7^ f2. Indeed 

£ a = - (d a $ + fid a #) u* + F ah u fc = (fi- ft)<9 a W +F a6 ti b = a, 6 =1,2; (19) 

thus, if F ah ^ 0, then Cl^Q. 
From eq. (|T6)) it also follows 

u 1 

*,2 = 2*,1» ( 2 °) 

which, differentiated with respect to ar, gives 



u 

*,i = -a5a*,i- ( 21 ) 



Using the continuity equation 

d 



u a a = -—HV^n) 1 (22) 



where n is the baryon number density, eq. (I21|) can be transformed as follows 



If we now define 



or \u z I j dr 



— ln(* 1) = u a 9 a ln(* = -w 2 ( - ) = -u% + — ^ = — In [nu 2 ^] . (23) 



we find 



* 1 

C^^^, (24) 

^-C = u a C a = Q => u 1 C 1 +u 2 C 2 = 0, (25) 
dr 

which, together with eq. (|^U]l implies that C is a function of $ only, i.e. C = C(^). 
By replacing $j = C^^nu 2 -^— g in eq. (fT9|) we find 

Fi2 = -(n-fi)Cnu*V=5. (26) 

Then, if $1 and do not coincide, F12 7^ 0; consequently, the magnetic field has both poloidal and toroidal components. 
A possible interpretation of f2 is the following (see for instance [l6|). From eq. (|26|) we find 

<> fi F|! 



from which we see that the fluid angular velocity Q has two contributions: the first, £l, can interpreted as due to 
the stellar rotation, the second is clearly due to the electromagnetic field. Although this interpretation is purely 
conventional, since we are considering a non rotating star, we shall assume Cl = 0, and consequently A Q = $ = 0. 
Thus the form of the vector potential is 

A^r,9) = {0,A r ,A 9 ,V) . (27) 
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3. Electromagnetic current and Lorentz Force 



The Lorentz force is denned as 

fn = F^r , (28) 

where the electromagnetic current J M is given by 

J" = ; -. (V^F^) . (29) 



The stress-energy tensor of a perfect fluid with an electromagnetic field is 

fluid ' em 

where 



T'" = T% id + T% (30) 



= (P + 1>K«" + P9" V . r e M m - ^ - \g" v F^F p ^ . (31) 

By projecting the equation T^f v — orthogonally to u M , we find the relativistic Euler equation in presence of a 
magnetic field: 

(p + p)a p + p >tl + u p u v p M - fn = , (32) 

where = u u u^ v . Let us now consider the <fi component of this equation. Under the stationarity and axisymmetry 
assumption it becomes 

(p + p)a tj) + u^u a p, a - = ; (33) 

being 

(Z0 = u^u^.fj, = u a U0, a - u^vyT^v = u a M0, Q + -u^u" = u a U0, a , (34) 
using the first law of thermodynamics, u a p^ a — £ i £ u a n Q , eq. (133(1 gives 



U = P -^u a (rax ) a . (35) 



If meridional currents are zero, — 0. 



B. The equations for the vector potential 

The background geometry of the star in coordinates (t, r, 8, <fi) is 

ds 2 = -e^dt 2 + e x{r) dr 2 + r 2 {d6 2 + sin 2 8d<j) 2 ) = gffldx^dx" (36) 
u ia ^ = (e-"/ 2 , 0,0,0), (37) 

where v (r), A(r) are the solution of Einstein's equations for an assigned equation of state. If a magnetic field is present, 
from the expression of F pu in terms of the electric and magnetic field — u^Ev — u u E p + e pva f3U a B° we see that, 
if = 0, then F^ = O(B); consequently, also the vector potential is of order O(B). By using a function A(r, 8) 
such that A t g = Ag, we can gauge away the 9 component of the vector potential (|2~7| . By introducing the function 

E(r, 8) = e~^~ (A r — A, r ), it then becomes 

A„ = (o,e^£,0,*) , (38) 

with £ and * of order O(B). 

The magnetic field induces motion in the fluid, and consequently induces a perturbation on the components of 
the four- velocity Su a , on the pressure and energy density (Sp and dp, respectively), and on the metric 5g pv . Since 
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^em = 0(i? 2 ), linearizing the equation T^ v = (and using the vanishing of the space components of u a when the 
magnetic field is absent), it is easy to see that Su a = 5p = 5p = 0(B 2 ). In a similar way, from the linearized Einstein 
equations it can be shown that Sg^v = 0(B 2 ). Thus, from eq. (|35|) we see that, since Su a — 0(B 2 ), the ^-component 
of the Lorentz force is fj, = 0(B 4 ) and for this reason hereafter we shall set it equal to zero. This condition will 
be used to further simplify the expression of the vector potential. We stress that the condition = comes from 
the fact that = 0(B 4 ), but we do not assume that meridional current are zero. If we compute from Maxwell's 
equations and impose = we find 

U = (E, ee + cot 0E, fl ) ¥, r - E,e r *, e = , (39) 

therefore, if we define 

f = sin0£,<>, (40) 

we have ^> g^ r — $ |t .$ g = 0; this equation implies — 4 r (^ r ), and consequently, since "J = O(B) and \& = 0(B), we 
can write \& = where £ is a constant of order 0(1). The equation 

C# = sm0£,0 (41) 

is satisfied by 

£ = (a 

* = sin 6a.fi, (42) 
with a — a(r,9). Thus, the vector potential can be written as 

= (0, Ce (A_l/)/2 a, 0, sin 6a <g ) . (43) 
As a consequence, the magnetic field takes the following form 

Bp = ^^(0,^ (smea.g) ,- (smdafi) r ,~Csm 2 de^/ 2 afi) . (44) 

From this expression we see that the coefficient £ (or the dimensionless quantity where R is the radius of the 
star) represents the ratio between the toroidal and the poloidal components of the magnetic field. Since, as discussed 
in the introduction, a magnetic field configuration with prevailing toroidal component is expected to be unstable, we 
will not consider configurations with £i? ^> 1. 

Assuming the form (|43[) of the vector potential, we find (neglecting the metric perturbations, which contribute to 
higher orders of B) 

f a = (wx9a, g ), a , J * , (45) 



r 2 sin 2 1 



where 



= - C 2 ^— sin 0a.fi . (46) 

47T 

We shall now show that f a can be written as (p + p) times the gradient of a function of (r, 8) . Let us consider the 
a-components of Euler's equation: 

(p + p)a a + p, a + u a u b pfi - f a = . (47) 

We remind that: 

u* = 6u* = 0{B 2 ), g 0i = 8g 0i = O(B 2 ), f a = F a ^ = 0(B 2 ) i = 1,2,3. 

Consequently, the term u a u b p t b is 0(B 4 ). We shall now compute a a and p jQ up to terms of order 0(B 2 ). 
The acceleration is 

a a = u^Ua-n = u b u a fi - uVr o(1 „ ~ ^u^g^a 

= 2 [("°) 2 foo,a + 2u u i g i,a + «*«%•,„] ~ i(w°) 2 5oo,a • (48) 
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We also find 

g^uPu" = -1 = (u°) 2 g 00 + 2u°u i g oi + uV 9ij ~ (u°) 2 g 0Q , (49) 

then 

a a = \ (M-9oo)), ■ (50) 
From the first principle of thermodynamics, written for a barotropic equation of state p = p(p), we find 

P,«=(p+J0(k— 1 ■ (51) 



If we introduce the function 



-P + P 



77 



X = In V^ra — , (52) 



using eqs. (|45)) , (|50|) and (fSTj) . eq. (|47|) becomes 



(p + p)x, a = (sin 0a,g), a „ j0 2/1 - (53) 



r 2 sin 2 1 



This equation is equivalent to eq. (12) of [l3|. From ([53| we find 

X,i2 ~ X,2i = (singo, 9 ),i ( , / ) - (sin fla, e ), 2 ( 2 . 2 lt | 7 ) = > ( 54 ) 

yr^sm 9(p+p) J 2 ^rsin 9(p + pj J 



hence 

J<4 



F(sin0a, e ). (55) 



r 2 sin 2 9(p + p) 
By expanding in powers of B we find 

-j^^iL—.^ + ^^ + otfl.,, (56) 

thus, the 4> component of the electromagnetic current can be written as follows 

J 4> = ( 2 ^— sin 9a fi + [cq + a sin 6a. e ] (p (0) + p(°) )r 2 sin 2 + 0(B 2 ) . (57) 

47T 

In the next section we will expand a{r, ff) in Legendre's polynomials; if we assume ci 7^ 0, different harmonic compo- 
nents of the field couple. Following [la [3 [3 Gil > hereafter we shall assume c\ = 0. With this simplification, 

J = ( 2 ^ sin 9a, e + c (p (0) + p m )r 2 sin 2 9, (58) 
where Co is a constant of order 0(B). The r, components of the current are simply: 

= — i y (sin 9a y e),e, (sin 9a t6 ). r . (59) 

47rsm0 \ r z J 

The electromagnetic current is the sum of two parts: 

j„ = j% + 4 ( 6 °) 

with 

JP = (0,0,0,cor 2 sin 2 ^(p (0) +p (0) )) 

4 = (61) 

J? is the source of the poloidal field (which does not depend on £); is the source of the toroidal field (proportional 
to C) and it is parallel to the magnetic field. Note that 
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• since J* oc £? M , it follows that F^ v J tv = 0; 

• when JP — (i.e. when Co = 0), then / p = F^J" — 0; therefore in this case the magnetic field is force free. 

If outside the star we assume there is vacuum, currents must vanish. As the poloidal current is proportional to 
p(°) -|- , it automatically vanishes; conversely, the toroidal current vanishes only if £ = 0, i.e. if the toroidal field 
vanishes. Therefore, in vacuum only poloidal fields (with no current) are allowed. 

If outside the star there is a magnetosphere, the situation is different because currents can be present, and conse- 
quently toroidal fields can exist. In any event, since the energy density in the magnetosphere is negligible with respect 
to that prevailing in the stellar interior, JP is negligible and the magnetic field is force free. 



C. The relativistic Grad-Shafranov equation 

If we expand the function a(r, 9) in Legendre's polynomials: 

oo 

a(r,e) =Y j ai(r)P l {6) , (62) 
(=1 

the vector potential (14"3")1 and the magnetic field (|4"4"| become 



i 



1 . n ( — A V r — \ r _ x 1(1+1) \ 

J<j, = -—sing y j Pifi I e aj >rr + - L — — -e aj, r ^ — a U ■ ( 65 ) 



A/2 

Bfj, = \0,-^l(l + l)aiPi.-e- x/2 ai, r Pi,e,-Ce~ u/2 aismePi,e) ■ (64) 

From Maxwell's equations we find = j^F^ , which gives 

V.r ~ A, r _ A /(/ + 1) 

Using the expansion ([62]) . eq. (l58|) gives 

J = C 2 — V a ; sin 6>Pi, e + c (/9 + p)r 2 sin 2 (9 = C 2 — > sin flails - c (p + p)r 2 sin <9P M . (66) 

4"7T 47T 

I I 

Notice that the poloidal current (i.e. the term in Cq) introduces an I — 1 dipole component. The linearized relativistic 
Grad-Shafranov equation [l6[ is found by equating eqs. (|B5]) and (frJo) (see also (l7j): 

e- A < + ^^e- A ai + f - ^ d = A <P + , (67) 



2 

e- A aJ + \X 2 e~ v ~ " v "^ ^ j a; = (I > 1) . (68) 

Hereafter we will consider only the solution of equation (f6"T|) corresponding to / = 1 , in which case the vector potential 
and the magnetic field (j6"4"|) become 

A M = (0,Ce (A_,y)/2 aiCos6l,0,-aisin 2 6i) (69) 

/ e A/2 \ 

Bf, = (0,-2—aicos9,e~ x/2 a 1 sm9Xe~ v/2 a 1 sm 2 9j . (70) 

It is convenient to express the magnetic field in terms of the orthonormal tetrad components (i.e. those measured in 
a locally inertial frame) in the background metric (|36p . i.e. 

B( r ) = — ^-cos# (71) 

B {g) = ^sin9 (72) 

B W = C -sin 9. (73) 
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D. Boundary conditions and matching with the exterior 

Different choices of the boundary conditions and of the matching conditions of the interior and exterior solutions of 
the Grad-Shafranov equation, correspond to different physical configurations. We shall consider the following cases. 

• Magnetic field extending thr oug hout the star. This configuration has been studied in the literature in several 
papers (for instance, in (l6l . Il7l|): however, it conflicts with the common belief that the neutron star core is 
superconductor. Actually, if the superconductor is of type II, the magnetic field extends throughout the star, 
but it has a very complicated structure (it is "quantized" in flux tubes). Thus, the smooth magnetic field we 
consider in this paper is a rough representation of such configuration. 

If we impose a regular behavior at the origin (which implies a\(r ~ 0) = agr 2 + 0(r )), for each pair of assigned 
constants ao,co the solution ai(r) is unique. 

• Crustal fields. If matter in the core is a type I superconductor, the magnetic field is confined in the crust, i.e. 
within 

r c <r <R, (74) 

where r c is the inner boundary of the crust and R is the stellar radius. We choose r c = 0.9 R. By imposing a 
regular behavior near r c , i.e. a\(r > r c ) = ao(r — r c ) + 0((r — r c ) 2 ), for each pair of assigned constants ao, Co 
the solution ai(r) is unique. 

We shall assume that outside the star there is vacuum, currents vanish and ( — (see eq. (J5HJ)). Equation (|67p then 
reduces to 

2M\ „ 2 M ,2 
1-— J< + — a' 1 -^a 1 = Q; (75) 

its general solution (decaying at infinity) is a pure dipole 

2M\ 2M 2M 2 ~ 



f \ 3/i 2 

oi r = 7rr 

u ' 8M 3 



In | 1 

r ) r 



(76) 



where the constant y, is the magnetic dipole moment in geometrical units. The corresponding magnetic field has the 
form 

e A/2 



B^ = 0, -2—o-ax cos 0, e~ A/2 a' 1 sin 0,0). (77) 



On the surface of the star, the function <zi(r) solution of eq. (|6"7| has to be matched with the exterior solution (|76p . 
imposing the continuity of a\ and a^. The ratio ao/co is fixed by matching the quantity a[/ai (which does not depend 
on n). Once this ratio has been determined, the constants ao, cq are rescaled by a common factor, which changes the 
constant /i (and then the global normalization of the field) by the same amount. We fix this constant by assuming 
that the magnetic field at the pole is B po i e — 10 15 G. In this way, for each assigned value of ( we determine ao and 
co- 
in previous papers on magnetized stars [HI , boundary conditions have been imposed in such a way that not only 
the toroidal, but also the poloidal component of the magnetic field vanishes outside the star; as a consequence, the 
parameter £ can take only a discrete set of values, a fact for which we do not see a reasonable, physical explanation. 

The matching conditions we impose at the boundaries are different, and should be considered as an attempt to 
better approximate realistic boundary conditions. Let us see why. We remind that outside the star we assume there 
is vacuum and ( — 0. By comparing ([TO")) . ([77)) we see that if we choose 01,0^ to be continuous across the stellar 
surface then B r , Bg are continuous. However, if ( 7^ inside the star, B^ is discontinuous because, having set £ = 
outside, it vanishes there. Such discontinuity corresponds to a surface current 

J^rf = U 0j -C-^ ax sin05(r - R), Oj . (78) 

A true neutron star is surrounded by a magnetosphere, where fluid energy density and pressure are small, but currents 
do not vanish. There, the magnetic field has both poloidal and toroidal components, and both match continuously 
across the stellar surface with their interior correspondent. The values of Oi, a'x which would ensure the continuity of 
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both components, would be different from those we choose by imposing — outside the star; however, with our 
choice at least we allow the poloidal field, which extends all over the space and decays as r~ l ~ 2 , to be continuous, 
whereas outside the star we switch off the toroidal component which extends only in the magnetosphere and tends to 
zero smoothly at its edges. A more precise characterization of the field behavior at the stellar surface would require 
the modelling of the magnetosphere, which is beyond the scope of this paper. 

The function ai(r), solution of eq. (|67|) . which describes the vector potential inside the star, can have zeros at some 
points r = f$; conversely, the exterior, vacuum solution (|76p never vanishes. This happens also in the Newtonian 
limit; for example, in [23j it has been shown that the solution of the equation corresponding to eq. (|67p with Co = 
is a linear combination of the spherical Bessel functions 

. , . sinx cosx 



cos a; sin x 

ni{x) = S"" 

(79) 

where x = £r. Such combination vanishes at given values of x. In this case £, which has the dimensions of an inverse 
length, can be interpreted as a sort of wavenumber of the solution. 

In our case also, though we set cq ^ 0, the location of the points where a\(r) vanishes depend on £. If r = f < R 
is a zero of a\, then B r (f) = (see eq. ([TOT ), and the magnetic flux is confined within the spherical surface r = f. 
This means that the field lines inside the star are defined in disjoint domains. Although we do not have a physical 
interpretation for this configuration, in our study we will not exclude this possibility. 



E. The ellipticity of the star 

The stellar deformation, which we determine by solving the perturbed Einstein equations given in Appendix [XI can 
be expressed in terms of the stellar ellipticity. In the current literature there are two different definitions of ellipticity, 
which correspond to two conceptually different quantities. The surface ellipticity, e sur f , is (TtI [l8l . [24| 

(equatorial radius) -(polar radius) 

e S urf = -, ; r — \ • (80) 

( polar radius) 

It describes the geometrical shape of the star. It should be mentioned that a slightly different definition has been 
employed in [25|, [26|,[27|, i.e. e sur f = J (e stir /) 2 + 2e sur /. The surface ellipticity describes the external appearance of 
the star. 

A different quantity is the quadrupole ellipticity, eg, which is a measure of the mass quadrupole of the star 0, [HI, [lj| : 

e Q = -f (81) 

where I is the mean value of the moment of inertia of the star iy, and Q is its mass-energy quadrupole moment. For 
a stationary, axisymmetric compact object, Q can be extracted by the far field limit of the metric [2||[28j]. Indeed, 
it is the coefficient of the l/r 3 P2(cos0) term in the expansion of c/oo m powers of 1/r and in Legendre polynomials 
Pi{0): 

.goo^ 2Q-^P 2 (cos#). (82) 

As discussed in [H, [29| , in the weak field limit the mass-energy quadrupole moment reduces to 

Q= [ P (r,6)r 2 P 2 {cos9)dV , (83) 
Jv 

where V is the star volume. In this limit, the quadrupole tensor of the axially symmetric star can be expressed in 
terms of Q: Qij = diag(— Q/3, — Q/3, 2/3 Q), and the quadrupole ellipticity can also be written in terms of the inertia 
tensor 
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In the general case, the quadrupole ellipticity is a measure of the entire stellar bulk deformation. 

Since eg and e sur f are quantities with different physical meaning, they are in general different. They coincide only 
in the case of a constant density star, in the Newtonian limit, as shown in Chapter 16 of 30]. 

It is worth stressing that the quadrupole ellipticity is the quantity that should be used to evaluate the gravitational 
emission of a rotating star; moreover, it has been used to study the spin-flip mechanism proposed by Jones and Cutler 



III. RESULTS 



In this section we present the results of the numerical integration of eqs. (|67|) . (|A6I) and (|A5|) . 

As a test, we have first ran our codes for the polytropic star used in [TJ|, endowed with mixed (poloidal and toroidal) 
magnetic field, which vanishes outside the star. Thus, we impose a\ — on the stellar surface r = R, and solve the 
eigenvalue problem to find the set of values Q for which this condition is satisfied. We have reproduced the values 
of Q given in Table I of [16[ for different values of the stellar compactness, with an accuracy better than 1%. The 
corresponding magnetic field profiles and stellar deformations (surface ellipticity and mass-energy quadrupole) are 
also in full agreement with (l6j |. 

Furthermore, we have integrated the equations for the models considered in Ref. [14j ; there, non rotating, magnetized 
stars with only poloidal fields have been modeled by solving numerically the full set of non-linear Einstein's equations; 
magnetic fields are either defined throughout the star, or confined in the crust. Following we introduce the 
magnetic distortion factor f3, given by 



M 2 



?T72 , (85) 



where M. is the magnetic dipole moment, related to magnetic field at the pole, B po ^ by 



47T 

M = B pole R 3 — . (86) 

Here, /io is the magnetic permeability. The normalization factor Mq is given by 

4tt GI 2 

M = w ■ 87 

Mo R 

With this normalization, the coefficient (3 is dimensionless. Moreover, as eg = 0(B 2 ) and M. = 0(B), (3 is nearly 
independent of B, and indicates to what extent a star can be deformed by the magnetic field. It is worth mentioning 
that the magnetic dipole M. defined in (f8"6"|) differs from the quantity /i defined in eq. (fTB")) . since eq. (|56")l has been 
derived in the context of Newtonian theory. As in 14j, we use the equation of state (EOS) of Wiringa Fix and 
Fabrocini [3l[ (WFF), and consider an M = 1.4 M© star. 

When the magnetic field extends throughout the star, we find (3 = 0.505, while the authors of [TJ[ find j3 = 1.01; 
when the field is confined to the crust, we find f3 ~ 5, while the authors of [3| find a very large value: j3 ~ 10 2 . 
However, if we compute (3 from the same equation, but using e sur f instead of eg, we find (3 = 1.01 when the magnetic 
field extends throughout the star, and (3 ~ 10 2 in the case of crustal fields, in agreement with [Til ]. 

The reason why, when crustal fields are present, the factor (3 computed using e sur f is much larger than that 
computed using eg, is the following. The crust contains a very small fraction of stellar matter therefore, although its 
deformation is large (because the field lines are squeezed in a small region), it does not induce a big change in the 
distribution of matter in the stellar bulk. As a consequence, e sur f ^ eg. 



A. Deformations induced by different magnetic field configurations 

We shall now study how the stellar deformations induced by a mixed (poloidal and toroidal) magnetic field depend 
on the field configuration. To describe matter in the stellar core we use the equation of state of Akmal, Pandharipandc 
and Ravenhall [32| (denoted as APR2); we choose a star with mass M = 1.4 M© and a radius R — 11.58 km. The 
magnetic field is normalized assuming that its value at the pole is B po i e = 10 15 G. For the different configurations 
discussed in Section III Dl we find the magnetic field structure, the surface and quadrupole ellipticities e sur f, eg, 
and the maximal values of the internal poloidal and toroidal fields, B 7 p nax and B T t nax . The equations for the stellar 
deformation and the procedure to compute e sur f , eg are described in Appendix [XJ 
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We stress that it is important to determine if the magnetic star has a an oblate or prolate shape, i.e. to determine 
the sign of cq; indeed, as suggested by Jones and Cutler @, 0], if e-Q < the star could change its rotation axis 
due to viscous forces ("spin flip") becoming an orthogonal rotator (with magnetic axis orthogonal to the rotation 
axis), and the process could be associated to a large gravitational wave emission. In this respect, it is also important 
to determine the absolute values of the allowed quadrupole ellipticities, because if the star rotates around an axis 
different from the magnetic field symmetry axis, it emits gravitational waves with amplitude (T3 | 

AC 

ho ^J±n 2 i\e Q \, 



and frequency vgw = ^/(27r), where is the angular velocity. 

It is worth stressing that the current upper bound on neutron star ellipticity, i.e. |e q| S, 10 -6 , is obtained by 
evaluating the maximal strain that the crust of an old and cold neutron star can sustain [33, [Hj]. However, a large 
deformation may be induced by the effect of strong magnetic fields in the very early phases of the stellar life, when 
the crust has not formed yet. These deformation may persist as the star cools down, leading to final configurations 
having an ellipticity larger than the above limit. Let us now discuss the two field configurations described in section 
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FIG. 1: The profiles of B( r ) evaluated at = 0, and of and evaluated at 8 — tt/2, are plotted as functions of the 

normalized radius inside the star. The magnetic field is defined through the whole star. Each panel corresponds to a value of 
C: C = 0.15 km" 1 in panel (a), C = 0.37 km" 1 in (b), £ = 0.40 km" 1 in (c) and ( = 0.30 km" 1 in (d). Panel (d) refers to a 
value of £ exterior to the ranges (|90[) : thus in this case the magnetic field lines are defined in disjoint domains (see text). 
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FIG. 2: The projection of the field lines in the meridional plane is shown for £ = 0.15 km -1 (upper panel, left), £ = 0.37 km -1 
(upper panel, right), £ = 0.40 km -1 (lower panel, left) and £ = 0.30 km -1 (lower panel, right). The dashed circle in the lower 
panel on the right, separates two disjoint domains. The magnetic field extends throughout the star. 

1. Magnetic field defined throughout the star 

If the magnetic field is non vanishing through the whole star (see Section III D[) . we find that a%(r) has no nodes for 
r < R in two ranges 

range 1 : (0 < C < 0.2915) (89) 
range 2 : (0.369 < C < 0.46) . 

If £ lies outside these ranges, the field lines inside the star are defined in disjoint domains, as discussed in Section Hi Dl 
In Figure Q] we plot the magnetic field components versus the radial distance, for r < R. B( r ) is evaluated at 6 = 0, 
.B(g) and Br^ are evaluated at 8 = tt/2. The plots are shown for four values of £ (in km" 1 ): £ = 0.15 in Fig. QJi) 
(this value is in the range 1), ( = 0.37 in Fig. [lb), £ = 0.4 in Fig. [lb), (both values are in the range 2) and C = 0.3 in 
Fig. HH), which is outside the ranges 1 and 2. Different values of £ give qualitatively similar behaviours. We see that 
when £ approaches the lower bound of range 2, as in Fig. [Tfc>), the field components become much larger than in the 
other cases. Of course they cannot be arbitrarily large, since they must not exceed the virial theorem limit B < 10 18 
G [l3[. This is a peculiar behaviour, which is not observed if one approaches the other bounds of range 1 and 2, 
either from inside or from outside. The reason for such behaviour is that the configuration with £ = 0.369 km -1 is 
a singular limit. It corresponds to a configuration in which ai(R) = 0, i.e. the magnetic field is confined inside the 
star and vanishes outside. This is inconsistent with the boundary condition we impose, i.e. B po i e = 10 15 G. Thus, 
this singular value is unacceptable. However, values of £ approaching this limit can be accepted, provided the virial 
limit is not violated. We mention that, as long as B is smaller than the virial limit, the stress-energy tensor of the 
electromagnetic field is smaller than that of the fluid, and the perturbative approach we use is appropriate. 

The field profiles shown in Fig. [Tfcl) refers to a case in which inside the star the field lines are defined in disjoint 
domains: indeed, they cannot cross the sphere r = 0.37 R since B r (r — 0.37 R) = 0. 

The projection of the field lines in the meridional plane is shown in Figure [5J the four panels refer to the same 
values of C considered in Figure [T] Fig.[2ji) corresponds to ( = 0.3, i.e. to the case of disjoint domains: field lines do 
not cross the dashed circle in the picture. 
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The ellipticities e sur f and eQ are plotted in Figure [3] as functions of £ € [0,0.5]. Continuous lines correspond to 
values of C inside the ranges ([50)) (no nodes inside the star), while the dashed lines correspond to values of £ for which 
there is a node inside the star. 

For small values of £ (i.e. if the poloidal field prevails) the star is oblate (e sur f.Q > 0). As ( increases, the toroidal 
part becomes more important and the star becomes prolate (e sur f t Q < 0). In other words, the toroidal field tends to 
make the star prolate, while the poloidal field tends to make it oblate; this behavior has already been discussed in 
the literature, see for instance [Tg|. For larger values of C the behavior is different. If we exclude values close to the 
singular point £ ~ 0.369 km -1 , we find that the ellipticity is 



w/,qI - io~ 6 - 10- 



(90) 



If we approach the value C = 0.369 from either sides, then e s 
virial theorem constraint, B < 10 18 G, corresponds to 



r f t Q < 0, and the deformation can be much larger; the 



-'surf 



< 2 x 10" 



|eq|<10 



-3 



(91) 



Thus, for a large range of values of £, the magnetic field induces a shape, either prolate or oblate, with | cq | ~ 
10~ 6 — 10~ 5 ; however, for very particular values of C, the star can have a strongly prolate shape (eQ < 0), with |eq| 
as large as 10~ 3 . 

In Table Q] we give, for selected values of £ in the range [0,0.5] km -1 , the surface and quadrupole ellipticities, and 
the maximal values of the internal poloidal and toroidal fields. It is interesting to note that for values of £ < 0.1, 
e sur f ~ 2eg. We find a similar behaviour when the ellipticity is induced by rotation and no magnetic field is present. 
Indeed, by integrating the equations of stellar deformation to second order in the angular velocity as in [27j , we find 
that, for a large variety of neutron stars EOS, e sur f ~ 2eQ for Q < 0.1 il ms , where Cl ms = \J M / R 3 . 



3e-05 



2e-05 




FIG. 3: Surface and quadrupole ellipticities as functions of £ for a star with mass M = 1.4 M©, and equation of state APR2. 
The magnetic field extends throughout the star. The dashed (solid) lines correspond to models for which ai(r) has nodes (has 
no nodes) inside the star. 



2. Crustal fields 



When the magnetic field is confined to the crust, we find that a\{r) has no nodes inside the star for 

< C < 1-085. 



(92) 



ai(R) for all values of £, therefore crustal field do not exhibit the singular behavior discussed in Section fill A II 

In the left panel of Figure [4] we show, for r < R, the profiles of Bf r ) evaluated at 6 = 0, and of Brg\ and 
evaluated at 6 = 7r/2), for £ = 0.5 km -1 . Different values of ( correspond to qualitatively similar behaviours. We see 
that the interior field is one order of magnitude larger than the surface field; this behavior, peculiar of crustal fields, 
is common to all values of £. The projection of the field lines in the meridional plane is shown in the right panel. 
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0.15 
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1.447 x lO" 6 


3.998 


2.948 


0.2 


-2.661 x 10" 7 
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0.25 


-2.643 x 10" 6 
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2.219 


3.303 


0.30 


-3.433 x 10" 6 


-5.343 x 10~ 7 


1.610 


3.264 


0.37 


-1.106 x 10" 3 
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518.0 


557.0 


0.35 


6.375 x 10~ 6 


4.273 x 10" 6 


21.52 


22.15 


0.40 


-2.220 x 10" 5 


-9.313 x 10" 6 


26.28 


29.05 


0.45 


-1.062 x 10" 5 


-1.263 x 10" 6 


18.89 


20.51 


0.50 


2.773 x 10 -6 


5.410 x 10 -6 


26.68 


27.50 



TABLE I: Surface and quadrupole ellipticities, and maximal values of the internal poloidal and toroidal magnetic fields, are 
tabulated for different values of 
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FIG. 4: The profiles of B( r ), evaluated at 6 = 0, and of and B^ evaluated at 8 — ir/2, are plotted for £ = 0.5 km 1 in the 
crust (left panel). The projection of the field lines in the meridional plane is shown in the right panel, for the same value of 



In Figure [5] we show the ellipticities as functions of £; continuous lines correspond to values of £ inside the range 
(f9"2)l . dashed lines to values outside that range. We see that, as discussed in the previous section, the geometrical shape 
of the star is oblate for small values of £ (for which e sur f > 0) and prolate for larger values: the surface ellipticity is 
a monotonically decreasing function of £. Conversely, the quadrupole ellipticity is always positive and, in modulus, 
much smaller than e sur f, even for values of £ larger than those considered in Figure [5] We note that these results rule 
out the Jones-Cutler mechanism in the case of crustal fields, since it can only occur when eQ < 0. As explained in 
Section IITEI the reason why e sur f 3> eQ is that, though the crust deformation is large since the field lines are squeezed 
in a small region, it does not induce a big change in the distribution of matter in the stellar bulk. 

In Table ILT1 we give, for selected values of £ in the range [0, 1.5] km -1 , the surface and quadrupole ellipticities, and 
the maximal values of the internal poloidal and toroidal fields. Comparing Tables HI and ILTl we see that, for crustal 
fields, typical values of e sur f are two orders of magnitude larger than for fields extending through the whole star. The 
quadrupole ellipticity is, typically, one order of magnitude larger: 

|e sur/ | - 1(T 4 - 1(T 3 , e Q ~ 1CT 5 - 1(T 4 , (93) 

with the exception of the models with £ close to 0.369, for which the deformation is larger in the case of fields extending 
throughout the star. 
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FIG. 5: The surface (left panel) and quadrupole (right panel) ellipticities are plotted as functions of £ for a star with mass 
M = 1.4 Mq and equation of state APR2, when the magnetic field is confined to the crust. The dashed (solid) lines correspond 
to models for which ai(r) has nodes (has no nodes) inside the star. 
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7.483 x 10~ 4 
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27.80 





0.2 
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0.4 
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23.29 
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0.6 


3.744 x 10~ 4 


1.310 x 10~ 5 


18.00 


5.616 


0.8 


1.313 x 10~ 4 


6.352 x 10~ 6 


11.19 


7.361 


1.0 


-1.259 x 10" 4 


3.113 x 10~ 6 


6.580 


9.130 


1.2 


-3.653 x 10~ 4 


4.798 x 10~ 6 


8.216 


10.91 


1.4 


-5.564 x 10~ 4 


1.127 x 10~ 5 


12.30 


12.70 


1.5 


-6.257 x 10" 4 


1.581 x 10" 5 


15.79 


13.59 



TABLE II: Surface and quadrupole ellipticities and maximal values of the internal poloidal and toroidal magnetic fields, are 
given for different values of (, in the case of crustal fields. 



B. Comparison between magnetic and rotational deformations 

Both rotation and magnetic field contribute to the ellipticity of the star, i.e. e sur f y Q — e^ ur j q + ef ur jq. It is 
interesting to compare the two contributions, evaluated in the range of parameters typical of observed magnetars 
(SGR and AXP), i.e. @ 

0.6 x 10 14 G < B < 7.8 x 10 14 G (94) 
5-2s < T < 11.8s, (95) 

where T is the rotational period. It should be mentioned that, as explained in |7j, only eg contributes to the spin-flip 
process, which occurs when eg < 0. 

We have computed e^ r jQ for an M = 1.4 M star with equation of state APR2, using the codes, developed by 
some of us [13], which describe the structure of a non magnetized, rotating star, up to 0(Q 3 ); ef ur tQ have been 
computed using the approach described in this paper. 

In Figure[H]we show |e sur /| and |cq | as functions of for the two magnetic field configurations described in Section 
III Dt field throughout the star (upper panels), and crustal fields (lower panels). The two solid lines correspond the 
|e^ tr ^g| computed for B po i e equal to the minimum and maximum values of the range ([M]) . The shadowed region 
corresponds to the rotation contribution, e^ ur ^ q, for rotation periods in the range (|95[) . The dashed lines correspond 
to T = 1 s and T = 0.1 s, outside that range and smaller than the observed periods of SGR's and AXP's: we show 
these values since they may possibly occur in young magnetars. From Figure [H we see that for the observed magnetars 
\ e Q surf I ^ s typically larger than |eg sur j\- This behaviour is magnified when crustal field are present (lower panels in 
Figure [S]). The rotational contribution may significantly exceed that of magnetic field only for stars rotating faster 
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(dashed lines). The solid line minima in the pictures correspond to the points where eg sur j = 0; there the effects of 
the poloidal and toroidal fields balance and the ellipticity changes sign. 
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FIG. 6: |e sur /| and |eg| are plotted as functions of (, for two magnetic field configurations: field throughout the star (upper 
panels), and crustal fields (lower panels). The two solid lines refer to |e^ r y q| computed for B po i e — B m in,max corresponding 
to the extrema of the range (|94[) . The shadowed region corresponds to e^, r <g evaluated for rotation periods in the range (|95fl . 
The dashed lines correspond to smaller rotation periods. 



C. Deformation of magnetized stars with different masses and EOS 

The results discussed in previous sections where obtained for a star with mass M — 1.4 Mq and equation of state 
(EOS) APR2. We shall now see how the results depend on the EOS and on the stellar mass. To this purpose, as an 
example we shall consider three different EOS: 

• APR2 [32| . derived within the non-relativistic nuclear many-body theory, assuming that the star is made of 
ordinary nuclear matter; the maximum mass is M max = 2.202 Mq. 

• G240 [35l | , derived within the relativistic mean- field theory and allowing for the presence of hyperons in coexis- 
tence with ordinary nuclear matter; M max = 1.553 Mq. 

• QS, based on the MIT bag model 36] (with B = 95 McV/fm 3 , a s = 0.4, m s = 100 MeV), assuming that the 
star is a bare quark star, i.e. composed entirely of deconfined quark matter; M max = 1.445 Mq. 

G240 with hyperons is a very soft EOS, QS is very stiff (for a comparative discussion of these EOS see refs. p7l. [37|). 
Furthermore we shall consider the two magnetic field configurations discussed in Section [II Dl and two values of mass, 
M = 1.20 Mq and M = 1.40 Mq. 
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In Figure[7]we show e sur f and cq as functions of £, for M = 1.2M® (upper panels) and for M = 1.4M Q (lower 
panels), for the selected EOS, when the magnetic field extends throughout the star. We see that, as expected, softer 
EOS and smaller mass correspond to larger deformations. For all masses and EOS, we find the same qualitative 
behaviour shown in Figure[3J 

In Figure [51 e sur f and eQ are shown in the case of crustal fields. We find that e SU rf depends strongly on the mass 
and on the EOS. For an assigned EOS, changing the mass from 1.4 Mq to 1.2 Mq, e sur f increases by a factor ~ 10, 
and eQ by factor ~ 2. If we fix the mass and change the EOS we find 

e surL Q{G240)/e surf , Q (APR2) ~ [2 - 4] , 

whereas 



e surf>Q (APR2)/e surf>Q (QS) ~ [10-100] . 

Finally, we find that for the stiffest EOS we consider (QS), for some values of £ the quadrupole ellipticity can become 
negative, whereas this never occurs for APR2 and G240. 

It should be stressed that when magnetic fields extend throughout the star the dependence of the ellipticities on 
the EOS and on the mass is considerably weaker (Figure [7]). 
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FIG. 7: Surface and quadrupole ellipticities as functions of for different equations of state, and magnetic fields extending 
throughout the star. The stellar mass is M — 1.20 Mq (upper panels) and M = 1.40 Mq (lower panels). 



19 



CD 




o 

CD 



0.002 
0.0015 

0.001 
0.0005 


-0.0005 
-0.001 

-0.0015 
-0.002 





M=1.4 M 






APR2 — ■ — 




G240 — «— 




QS — — 





0.2 0.4 0.6 o.; 

C (km" 



1 1.2 1.4 




0.2 0.4 0.6 0.1 

C(km" 1 ) 



1 1.2 1.4 



0.0001 



5e-05 



-5e-05 



-0.0001 




0.2 0.4 0.6 0.8 1 1.2 1.4 

Ukm' 1 ) 



FIG. 8: Surface and quadrupole ellipticities are plotted as in Figure [7] in the case of crustal fields. 



D. Higher order multipoles 



In this paper we have focused on the study of dipole (I = 1) magnetic fields, which decay as r 
dominate far away from the star. In this Section we solve the Grad-Shafranov equation (|6"8)) for I = 
poloidal and toroidal components: 



— A // 
e a 2 



A' 



e a' 2 



6 



a 2 =0. 



For fields extending throughout the star, we impose a regular behaviour near the origin: a 2 (r ~ 0) = a^r c 
In the case of crustal fields the regularity condition is imposed near the crust-core interface, i.e. a 2 (r 
a 2 {r - r c ) + 0((r - r c ) 2 ). 

We assume that the field vanishes outside the star, i.e. a 2 (r > R) = 0. Continuity of (and then of B,, 



If 



' 2 and therefore 
2, including both 



(96) 

+ 0(r 5 ). 

> r B ) = 

.) on the 



), to select the discrete 
which corresponds to 



stellar surface, implies a 2 (R) — 0; thus we have to solve an eigenvalue problem (like in 
set of values £ = Q for which the boundary conditions are satisfied. The eigenfunction a 2 (r) 
£i, has i-nodes, one of which is located at the stellar surface. We note that, since as mentioned in Section III C[ only 
the current contributes to a; when I > 1, we do not have as much freedom as in the / = 1 case, when we used the 
constant cq to impose a' 2 {R) = 0; consequently, a' 2 is discontinuous (and so is Bg) on the stellar surface. 

For a star with M = 1.4 M Q , described by the EOS APR2, we have determined the field configurations corresponding 
to the first five eigenvalues Q. Since the field vanishes on the stellar surface, we normalize B by choosing a$ such that 
the maximum value of the magnetic field inside the star is B max = 10 16 G, i.e. of the same order of magnitude of the 
field considered in Section IIII Al for I = 1 (see Tables HHl]) . Then we have solved the equations of stellar deformation 
given in Appendix (|A 2[) . finding the surface and quadrupole ellipticities. 

The projection of the field lines in the meridional plane is shown in figure [5] for the first two eigenvalues. The upper 
panels refer to fields defined throughout the star, the lower panels to crustal fields. 
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FIG. 9: The projection of the I — 2 field lines in the meridional plane is shown for £ = £i (left panels), and f = £ 2 (right 
panels). The upper panels refer to fields extending throughout the star, the lower panels to crustal fields. 
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TABLE III: First eigenvalues and the corresponding surface and quadrupole ellipticities, for I — 2 magnetic fields extending 
throughout the star. 
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TABLE IV: First eigenvalues C,i, and the corresponding surface and quadrupole ellipticities, for I — 2 crustal fields. 



In Tables IIIII and IIVI we give the first five eigenvalues d and the corresponding ellipticities, for fields extending 
throughout the star and for crustal fields, respectively. In the first case the ellipticities are always negative and of the 
order 

\e SU rf, Q \ ~nr 7 - 10- 6 (97) 

i.e. smaller than for the I = 1 fields. For crustal fields, ellipticities are always positive and of the same order of 
magnitude as for I — 1, i.e. 

e S urf, Q ~ 10" 5 . (98) 
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IV. CONCLUDING REMARKS 



In this paper we solve Einstein-Maxwell's equations, using a perturbative approach, to study the structure of the 
magnetic field of magnetars, and to find the deformation it induces on the star. We extend previous works on the 
subject fl q . fl7l . jll . in several respects: we include toroidal fields inside the star, thus removing the assumption, 
used in [13lll4lll5|. of circular spacetime; we determine both the surface ellipticity and the quadrupole ellipticity; we 
explore various field configurations, corresponding to different boundary conditions; we compare the effects produced 
by a magnetic field and by rotation on the stellar structure; we study how different equations of state and masses 
affect the magnetic field structure and the quadrupole ellipticity it induces; we solve the equations for higher order 
(I = 2) multipoles. 

In summary, the main results of our study are the following. 

• Crustal fields induce surface deformations much larger than fields extending throughout the star, but the 
quadrupole deformations are comparable in the two cases. Typically, crustal fields produce oblate, rather 
than prolate shapes. 

• For particular values of the parameter £, representing the ratio between toroidal and poloidal components, the 
magnetic field inside the star and the deformation can be extremely large; such configurations correspond to 
prolate shapes. 

• For the typical rotation rates of observed magnetars, the deformation induced by rotation is much smaller than 
that induced by the magnetic field. 

• Neutron stars with the same magnetic field, but with softer EOS or smaller mass, exhibit larger deformations. 

• If the magnetic field extends throughout the star, the deformations induced by higher order (I = 2) multipoles 
are one order of magnitude smaller than the dipolar contributions; for crustal fields, they are comparable. 

As a future extension of this work, we plan to study the effect of couplings between different multipoles, which we have 
neglected in the present paper, and to determine their relative weights. Furthermore, the equilibrium configurations 
we have found will be used as background models to study the oscillations of highly magnetized neutron stars. 



Acknowledgments 

We thank Nils Andersson, Juan Antonio Miralles, Luciano Rezzolla, Lars Samuelsson, Kostas Glampedakis and 
Riccardo Ciolfi for useful suggestions and discussions. 



APPENDIX A: THE DEFORMATIONS OF THE STAR 

The metric of a non rotating star deformed by a magnetic field can be written, up to 0(B 2 ), as [l6[ 

ds 2 = -e p {l + 2[h + h 2 P 2 {cos6)])dt 2 + e x (^ + ^-[m a + m 2 P2(cos9)]^ dr 2 

+r 2 [l + 2k 2 P 2 (cos9)} (de 2 + sin 2 9d(t) 2 ) 

+2 [i^ (cos 0) + i 2 P 2 (cos6) + i 3 P 3 (cos0)}}dtdr 

(d d d "N 

vi qqPi (cos 6>) + v 2 — P 2 (cos 0) + v 3 — P 3 (cos 0) j dtdxj> 

( d d \ 

+2sin<9 l W2 -P 2{cos9 )+ W3 -P 3 (cose))drdt. (Al) 

The perturbed quantities {hi{r), rrii(r), TOj(r), fcj(r)) (i=0,2) and (ii(r), Vi(r) ii(r)) (i=l,2,3) are found by solving 
the linearized Einstein equations 

SG^ = 8nST^ . (A2) 
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The pressure p and the energy density p can be expanded as p = p^ ' + p = p(°) + <5p, with 

*p(r, 0) = (6p + 5 P2 P 2 (cos6)) 

p (0)' 

= j m j( 5 Po + Sp2P2(cos6)). 



(A3) 
(A4) 



1. Deformation induced by a dipole (/ = 1) magnetic field 



As discussed in section ITTC| a(r, 8) — a\(r)P\(cos8); by expanding the components (rr), (r9), (88) — sin 2 8(4>4>) of 

in spherical, tensor harmonics and by considering the I = 2 equations, which 



the perturbed Einstein equations (|A2 
give the stellar deformation, we have 
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Assuming regularity of h 2 and j/2 as r — > implies that near the origin 

fta =s ^ 2 , y 2 ~ Br 4 , 

where 



„ „ 16 9 
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^)-fa c (,(0) + ,(o)) + ^ 



(A5) 

(A6) 

(A7) 
(A8) 

(A9) 



It is worth mentioning that the terms in ai,a[ which appear in the definition of y 2 (eq. IA7[) are important: if they 
are not included (i.e. if we define y 2 = h 2 + k 2 ), the asymptotic behaviour (|A8[) is not satisfied. 
The quantities h 2l y 2 inside the star can be decomposed as follows: 

hi = cih 2 + h\ 

V2 = c iy %+y p 2 . (A10) 

For magnetic fields extending throughout the star, h 2 and y 2 can be found by integrating (|A5[) . (|A6|) from r — with, 
for instance, A = 1 and -B given by (|A9[) : ^2 an d 2/2 are the solutions of the associated homogeneous equations (i.e. 
with a\ = a'i = 0). 

When the magnetic field is confined to the crust, in the core (defined conventionally by < r < r c ) a\ = 0, and 
eqs. (|A5I) . (|A6|) are homogeneous; thus in this region h 2 = h 2 , y 2 — y 2 . On the crust-core interface r = r c , we impose 
ai = 0, a[ = const. We integrate the non-homogeneous equations starting at r = r c with the initial conditions 
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2V 
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2/2 M 
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>i(0) 2 . 



(All) 



The non- vanishing value for y 2 (r c ) follows from the requirement of continuity of h 2 + k 2 at the crust-core interface 
(see eq. (|A7| V 

The constant c\ in (j A10|) can be determined by matching the solution inside the star with the analytical solution 
in vacuum 1171: 



h 2 = KQ 2 2 (z) + h 2 (z) 
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(A12) 
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Here K is an integration constant, ai(r) is given by eq. (|76 
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QJ^ are the associated Legendre functions of the second kind 
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(A13) 

(A14) 
(A15) 

(A16) 

(A17) 



We have checked, both analytically and numerically, that (j A12[) is actually solution of (| A5|) . (|A6[) in vacuum. Matching 
hi and 1/2 at r = R allows to fix the constants c\ and K . 

The integration constant K is related to the mass-energy quadrupole moment of the star (see section III Ep by the 
relation 
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5Q 3^ 2 



Indeed, the asymptotic limit of h^ir) for r — > 00 is 
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Finally, we can compute the surface ellipticity of the star ([50)1 following the definitions of [l?], 124 



^■sur f 



3 /<5r 2 



= — -k 2 



2 V r 



3 / ^2 



2 V rp(°)' 



- = - 



2co&i 3/i2 3A)2 



(A18) 



(A19) 



(A20) 



r=-R 



where <5p = SpiPi, Sr — <5r;P; and 
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The relation (|A21|) is a consequence of Euler's equation. Indeed, from equations (|52[) . (f53"| . (f5fi)) . it follows that 

-P + P 



In V-.900- 



co sinfla^ + const. 



(A22) 



If we perturb (|A22[) . using the following relation which holds for a barotropic EOS 

Sp = nS 1 



we find (|A"2i"|) . 
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2. Deformations induced by a quadrupole (/ = 2) magnetic field 



We assume a(r,9) = a 2 (r)P 2 (cos9), and expand in spherical, tensor harmonics the (rr)-, (r6 
components of the perturbed Einstein equations (|A2|) . We find that the I = 2 equations are: 
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where we have defined 
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— a 2 a 2 - ^a 2 . 



Assuming regularity of h 2 and y 2 as r — > implies that near the origin 

h 2 ~ Ar 2 , y 2 - Sr 4 , 

where 

S = -2vrA f p(°) + ^ I • 



(A23) 



(A24) 



(A25) 



(A26) 



(A27) 



The integration of eqs. (|A23|) . (|A24|) and the determination of Q and e SU rf can be performed as in the previous 
section. 
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